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ABSTRACT. Mass-balance analysis of the Greenland ice sheet based on surface elevation changes 
observed by the European Remote-sensing Satellite (ERS) (1 992-2002) and Ice, Cloud and land Elevation 
Satellite (ICESat) (2003-07) indicates that the strongly increased mass loss at lower elevations 
(<2000 m) of the ice sheet, as observed during 2003-07, appears to induce interior ice thinning at 
higher elevations. In this paper, we perform a perturbation experiment with a three-dimensional 
anisotropic ice-flow model (AIF model) to investigate this upstream propagation. Observed thinning 
rates in the regions below 2000 m elevation are used as perturbation inputs. The model runs with 
perturbation for 10 years show that the extensive mass loss at the ice-sheet margins does in fact cause 
interior thinning on short timescales (i.e. decadal). The modeled pattern of thinning over the ice sheet 
agrees with the observations, which implies that the strong mass loss since the early 2000s at low 
elevations has had a dynamic impact on the entire ice sheet. The modeling results also suggest that even 
if the large mass loss at the margins stopped, the interior ice sheet would continue thinning for 300 years 
and would take thousands of years for full dynamic recovery. 


INTRODUCTION 

In contrast to a near-balance status in the 1990s (Zwally and 
others, 2005), the Greenland ice sheet (GIS) has experi- 
enced a marked mass loss (1 71 Gta -1 ) since the early 2000s 
(Zwally and others, 2011). Substantial losses occurred at 
low elevations due to extensive surface melt and dynamic 
thinning at the margins (Zwally and others, 201 1 ), particu- 
larly in the outlet glaciers as a result of accelerated flow 
(Howat and others, 2007; Pritchard and others, 2009) and 
surface melting (Mernild and others, 2009). Although the 
inland portion of the GIS continued to grow from the 1 990s 
due to the increased accumulation rate, the mass gain 
above 2000 m elevation has decreased since the 2000s 
(table 1 of Zwally and others, 2011). After correction for 
enhanced firn compaction due to surface warming (which 
lowers the surface elevation) and bedrock motion, the 
approximate change in thickness of the ice column (d//df) 
was nearly unchanged during 2003-07 compared with the 
1990s. However, an increase in surface elevation (d H a CA / 
df) at higher elevations due to the increase in accumulation 
was approximately balanced by an increase in dynamic 
thinning (d/T bd /df) that also extended to higher elevations 
(fig. 10 in Zwally and others, 2011). Those results suggest 
that increased melting and strong dynamic thinning 
initiated at the margins propagates rapidly into the high 
elevations of the ice-sheet interior in some drainage 
systems (Fig. la). 

Previous modeling studies have found that the thinning 
and acceleration initiated by the change in boundary 
conditions at the ice terminus propagate inland (Parizek 
and Alley, 2004; Payne and others, 2004; Price and others, 
2008a, b, 2011; Nick and others, 2009; Joughin and others, 
2010). Most studies focused on the locally fast flow region or 
the dynamic changes in a short time period. 

In this paper, we examine the time response of the GIS to 
elevation perturbations at the margins using a three-dimen- 
sional (3-D) ice-flow model that incorporates anisotropic ice 


flow. The perturbation experiment is conducted using rates of 
thinning at lower elevations derived from the European 
Remote-sensing Satellite (ERS) and Ice, Cloud and land 
Elevation Satellite (ICESat) as the perturbation input to 
investigate the impacts of the recent strong mass losses at 
the margins on the interior of the ice sheet. 

Numerical model 

A 3-D time-dependent Anisotropic Ice-Flow model (AIF 
model) incorporating anisotropic ice flow was developed, 
which fully couples dynamics and thermodynamics. This is a 
higher-order model with longitudinal and vertical shear 
stresses. The model calculates ice-sheet thickness, distribu- 
tions of stress, velocity, temperature and flow properties 
throughout. In Cartesian coordinates (x, y, z), i denotes x 
and y, and j denotes z. The ice thickness is calculated based 
on the continuity equation: 

d ^=-V-(V i H)+B-M, (1) 

where ^ is the change in ice thickness H with time t. B is 
the net surface mass balance taken as the present-day 
constant (updated from Zwally and Giovinetto, 2000). M is 
the basal melting rate. V,- is the depth-averaged horizontal 
velocity, V, : 


Vi = 2 [ ijj dz + V b i. 
J 0 

(2) 

ijj is the shear strain rate. The basal sliding velocity V b and 
basal melting rate M are calculated based on the basal shear 
stress, T b = -pgHa, as 
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X 

II 
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. , V b T b 
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pLf 

(4) 


where A s is a constant (Wang and others, 2002a), Z* is the 
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Fig. 1 . Greenland ice sheet thickness changes, (a) Dynamic thickness changes from ICESat observation (fig. 7b in Zwally and others, 201 1), 
dHt,d/df, include both ablation and dynamic terms below the equilibrium line and only dynamic terms above the equilibrium line. The 
stronger propagation regions where the dynamic thinning extends inland in drainage systems 3 and 7 are indicated, (b) Model perturbation 
input: the difference in dynamic thickness changes (<2000 m elevation) between ICESat (2003-07) and ERS (1992-2002) observations, 
(c) Modeled results after 1 0 year perturbation, showing agreement with the observations in (a) that enhanced inland thinning occurs mainly 
in the western, eastern and southeastern regions of the ice sheet. Results for four locations along the flowline (dashed line) marked a, b, c 
and d are presented in Figure 3. The 2000 m contour is shown as a dotted line in (a) and a thin solid line in (b) and (c). The equilibrium line is 
shown as a dashed line in (a) and a thick solid line in (b) and (c). 


height above buoyancy, L f is the latent heat of fusion of ice 
(Price and others, 2008b), p is the density of ice, g is the 
acceleration due to gravity and a is the surface slope. The 
parameters are presented in Table 1 . V b and M occur when 
the basal temperature reaches the pressure-melting point 
(Eqns (9) and (10)). 

A successful ice-sheet model relies on an accurate specifi- 
cation of the flow law that governs the flow behavior of ice. 
Most ice-sheet models use a form of Glen's flow law (Glen, 
1958), which is not appropriate when anisotropy of the ice 
crystal fabric develops in the deeper parts of the ice sheet. 
Using the flow law for isotropic ice (Glen's law) leads to a 
flow velocity three to eight times less than observed values. 
Most ice-sheet models simply use a constant enhancement 
factor of 3 or 5 (e.g. Huybrechts and others, 1 991 ; Hvidberg 
and others, 1997) to account for this discrepancy. Some 
attempts have been made to describe the anisotropic ice flow 
by incorporating the effects of anisotropy into the constitutive 
relation (e.g. Azuma and Goto-Azuma, 1996; Greve and 
others, 2009). Alternatively, Glen's flow law can be modified 
to include anisotropic effects by introducing an enhancement 
factor defined as the ratio of the strain rate for anisotropic ice 
to the strain rate for isotropic ice, where an increase in the 
enhancement factor is associated with a strengthening of a 
near-vertical single-maximum fabric (e.g. Dahl-Jensen, 1 985; 
Wang and Warner, 1999). 

To provide reliable estimates of evolving ice-sheet geom- 
etry, the treatment of the flow of ice needs to be based on a 


proper understanding of ice rheology. In this model, we use 
the flow law from a previous study (Wang and Warner, 1 999): 

^-f(A c )A(T)r,^- 1 +cr 1 ) = 0. (5) 

Enhancement factor E( A c ), derived from laboratory measure- 
ments of ice rheology (Li and others, 1 996), accounts for the 
effect of anisotropy for ice: 

f(A c ) = Es(Ec/Es)\ (6) 


Table 1 . Physical parameters used in the model 


Symbol 

Parameter 

Value 

Unit 

A) 

Temperature parameter (T< 263) 

1.14 x10“ 5 

Pa" 3 a" 3 

A> 

Temperature parameter (T> 263) 

5.47 x 10 10 

Pa -3 a -3 

A s 

Sliding parameter 

2.0 x 10 -11 

nr Pa -3 a -1 

c 

Specific heat of ice 

1 52.5 + 7.12 T 

1 kg-’ IC’ 

g 

Gravity acceleration 

9.81 

ms 2 

k 

Thermal conductivity of ice 

9.828e~ 00057r 

W n-f 1 K~’ 

L, 

Latent heat of fusion of ice 

333.5 

k) kg 1 

n 

Flow law exponent 

3 

- 

Q 

Activation energy (T < 263) 

60 

kj mol 1 

Q 

Activation energy (T> 263) 

139 

kjmof 1 

R 

Universal gas constant 

8.314 

| moP 1 IC 1 

P 

Density of ice 

910 

kg m -3 

G 

Basal geothermal heat flux 

-0.07 

Wnf 2 
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where 

A c = Til/ \Jrfj + tI or A c = ini yje\ + £?■ ( 7 ) 

is a function of stress configurations. T s =10 and E c = 3 are 
respective enhancement factors for shear or compression 
alone. As A c varies with ice-sheet depth from 1 to 0 as the 
stress situation varies from vertical compression (or longi- 
tudinal extension, i.e. r,y= 0) to simple shear (i.e. r,-, = 0), the 
enhancement factor increases from 3 to 1 0, corresponding to 
ice fabric changes from a small-circle girdle pattern to a 
single-maximum pattern. Previous applications of this 
enhancement factor have successfully described the flow 
behavior of ice (e.g. Wang and Warner, 1999; Wang and 
others, 2002a; Breuer and others, 2006; Price and others, 
2007). 

The temperature parameter A(T) is taken as the Arrhenius 
relation (Paterson, 1994): 

A( T) = A 0 exp Xj'j , (8) 

where A 0 is a constant, R is the universal gas constant and 
Q is the activation energy for ice creep (Table 1). By 
assuming that the temperature gradients in horizontal 
directions are very small compared with the vertical 
gradient, the temperature Tis calculated from the equation 
for heat transfer: 


dT _ k dT 2 1 dk fdT\ 2 
dt pc dXf + pc dT l^cAY , ) 


— * (y I 

VrVT ~ w Wj + q ' (9) 


where c is specific heat of ice and k is the thermal 
conductivity of ice (Table 1). q = 2r,y£,y/pc is the internal 
heating due to the deformation. The vertical velocity w can 
be obtained from the horizontal velocities (Eqn (2)) based on 
the assumption that ice is incompressible. 

Two boundary conditions are required to solve Eqn (9). At 
the surface, the temperature is taken from the model input. 
At the base, the temperature is determined as 


T — T pmp 
dT_ G 
dX 7 ~~k 


( T > T pmp ) 
(T < T pmp) 


( 10 ) 


T p mp is the pressure-melting point. C is the basal geothermal 
heat flux (Table 1). 

The shear stress r, ; is calculated in terms of ice density p, 
the acceleration due to gravity g, the depth h and the surface 
slope a,: 


Tij = —pghaj. 


( 11 ) 


In Eqn (5) the longitudinal strain rate (in) is involved in the 
shear strain rate (e,y) to take account of the variations of flow 
properties and stresses along ice flow. The longitudinal strain 
rate is the gradient of the horizontal velocity 


It can be expressed as 

\ i 

dX; ~ dX-, 


L ./0 


d 

dXi 


dVi 
£ii ~dX : ' 

2f(A c )A(T)r,y dz 

[ Z 2E(\ c )A(T)T? i -\ ij dz 

Jo 


( 12 ) 


(13) 


which indicates that the flow law (Eqn (5)) also includes the 
effect of the shear stress and the longitudinal stress gradients 


which are considered to play an important role at some 
places in the ice sheet, such as the areas of rough bedrock, 
fast ice flow or ice divides (Budd, 1971; Pattyn and others, 
2008). If in is not included in Eqn (5) the flow law simplifies 
to a shallow-ice approximation form (flutter, 1983): 

eij = E(X c )A(T)rlj. (14) 

In this study we used the anisotropic ice flow law (Eqn (5)) 
and Eqns (1-12) to solve the ice thickness, shear and 
longitudinal strain rates, enhancement factors, temperatures, 
velocities, shear and longitudinal stresses. 

PERTURBATION EXPERIMENT 

The perturbations are taken from the difference of the ice- 
thickness changes (d/Tbd/df) below 2000 m elevation 
between ICESat (2003-07) and ERS (1992-2002), which 
represent the thickness change since 2003 (Zwally and 
others, 201 1 ). d/T bd /df includes the combined dynamic and 
ablation terms (d/T bd /df = dH d /dt+ dH b /dt). In the ablation 
zone, d/Tbd/df includes contributions from both ice dynam- 
ics (d/Td/df ) and melt (dHb/dt). In the accumulation zone, 
d/Tbd/df contains only the dynamic term (Li and Zwally, 
201 1 ). The decrease of d/Tbd/df above 2000 m is considered 
to be a result of the inland propagation from the mass loss at 
lower elevations of the ice sheet. At some locations (mainly 
in the north) the perturbations are small and positive, which 
would cause dynamic thickening. To emphasize the stronger 
dominant impacts of the mass losses, we only use the 
gridpoints with negative perturbations. All positive perturba- 
tions are set to zero (Fig. 1b). 

The experiment is intended to answer two questions: 
(1) How far and how rapidly will the strong mass loss from 
the coastal regions propagate inland? (2) How many years 
will the impact remain? 

We set up the modeling runs in three steps. The first step 
is to establish a steady-state ice sheet which corresponds to 
present-day conditions (Fig. 2), using the present-day basal 
topography (Bamber and others, 2001), surface mass 
balance (updated from Zwally and Giovinetto, 2000), 
surface temperature and balance velocity. Derived balance 
velocity is based on the ICESat observed thickness (http:// 
nsidc.org/data/nsidc-0304.html) and present-day surface 
mass balance using the scheme from Budd and Warner 
(1996). To establish the present-day steady-state ice sheet 
(Fig. 2b and d), the model iterates on the governing 
equations (Wang and Warner, 1999), with the balance 
velocity (Fig. 2c) as a target. To compensate for any 
overestimated velocity without tuning any parameter in the 
flow law, we adopted the scheme from Wang and Warner 
(1 999) by discounting integration of shear strain rate near the 
base of the ice sheet to match the balance velocity (Fig. 2c). 
This treatment produces strain-rate profiles with a band of 
high shear above the bedrock which can be compared with 
observations (e.g. Wang and Warner, 1 999, fig. 3; Wang and 
others, 2002b, fig. 4). 

In the second step, we put in the perturbations (Fig. 1b) 
and run the model for 10 years. The modeled results are 
compared with the observations (d/T bd /df). The third step is 
set up to answer the question of how many years the impact 
is sustained after the perturbations are switched off by 
running the model forward lOOka. The surface mass 
balance and temperature remain unchanged, which allows 
us to examine the impact from an isolated 10 year 
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Fig. 2. I nitial state of the Greenland ice sheet for perturbation experiment, (a) Thickness derived from ICESat observation (http://nsidc.org/ 
data/nsidc-0304.html). (b) Modeled thickness, (c) Derived balance velocity based on the ICESat observed thickness and present-day surface 
mass balance, (d) Modeled depth-averaged horizontal velocity. The 2000 m contour and equilibrium line are shown as a dotted line and 
dashed line, respectively. 
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perturbation | propagation 



Surface elevation (m) 

Fig. 3. (a) The rate of ice-thickness change averaged in 400 m 
elevation bands over the ice sheet versus elevations from the 
observations and the modeled results after 10 year perturbation 
(Fig. 1c) over the whole ice sheet (CIS, black solid line) and 
drainage system 3 (DS3, red dashed line) which shows the 
stronger propagation, (b) Corresponding depth-averaged horizontal 
velocities. The ranges of perturbation (<2000 m elevation) and 
propagation (>2000 m elevation) are shown by the arrows. 


perturbation of the thickness. The model does not include an 
ice-front calving process. dHbd/df is small and the total 
thinning for a 10 year perturbation at any perturbation 
gridpoint is a small portion of the thickness. Therefore, we 
treat the marginal ice edge as steady-state during the 
experimental period. 

The model has a horizontal resolution of 50 km, 
compatible with the resolution of the observations, and 31 
evenly spaced vertical layers. The time-steps are 1 year. The 
physical parameters used in the model are listed in Table 1. 
A finite-difference technique is used to solve the equations. 

RESULTS AND DISCUSSION 

The modeled results after 10 years of perturbation input 
show that the observed strong mass loss at low elevations 
certainly causes interior ice-sheet thinning (Figs 1c and 3a). 
The dynamic propagation penetrates deep into the ice sheet, 
to almost the highest elevations. The modeled thinning 
pattern (Fig. 1c) over the ice sheet is in agreement with the 
observations (Fig. la), showing that enhanced thinning 
mainly occurs in the west and east, as well as in the 
southeast. In contrast, the thinning is least in the north and 
southwest. This implies that satellite-observed thinning (or 
reduced thickening) in the higher elevations of the ice sheet 
is likely induced by the extensive mass loss at the margins 
during the past decade or so. 

The stronger propagation (Fig. 3a) corresponding to the 
higher horizontal velocities (Fig. 3b) in drainage system 3 
indicates that faster ice flow enhances the propagation. The 
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Fig. 4. Modeled ice-thickness changes with time (log scale) at four 
selected locations along the flowline in West Greenland (see 
Fig. 1c): (a) at the perturbation region; (b-d) at ~100km (b), 
~200 km (c) and ~400 km (d) upstream from the perturbation 
region. The insets in (c) and (d) plot the time using a linear scale to 
show the time when the maximum thinning rates are reached. 


lack of resolution of basal troughs at 50 km grid spacing may 
result in reducing the speed of the propagation, as the deep 
troughs lead to higher driving stresses and warmer ice and 
therefore faster ice flow and stronger propagation. Improve- 
ment of the modeling spatial and temporal resolutions may 
enhance the propagation and reduce the discrepancies in 
the thinning rates between the modeled results and the 
observations (Fig. 3a). 

After the perturbations are switched off, the ice sheet in 
the perturbation region starts to grow at a very high rate and 
then slowly returns to the previous steady state (Fig. 4a), 
while the ice sheet in the propagation region keeps thinning 
before recovering (Fig. 4b-d). The maximum thinning rates 
occur after 1, 55 and 350 years (i.e. years 11, 65 and 360 in 
Fig. 4b— d) at 100, 200 and 400 km upstream, respectively, 
showing a propagation wave: the farther from the perturb- 
ation region, the later the response occurs. The thickening 
rates during the recovery period are low compared with the 
ice loss rates during the period of thinning, so the ice sheet 
takes longer to return to its pre-perturbation state than it did 
to reach the thinnest state. 

The changes in ice-sheet mass through time (Fig. 5) show 
the impact of the 10 year perturbations. In the first 10 years, 
the interior of the ice sheet (propagation region shown as red 
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Fig. 5. Modeled ice-sheet mass change over the whole ice sheet 
(black line) and the ice sheet above 2000 m elevation (red line) as a 
function of time (a) in log scale for 1 0 ka and (b) in linear scale for the 
first 400 years, showing the times when the maximum mass losses are 
reached. The left-hand labels indicate the mass change relative to the 
maximum mass loss. The right-hand labels indicate the total mass 
loss (Ct) from the steady state for the whole ice sheet (black numbers) 
and for the ice sheet above 2000 m elevation (red numbers). 


line) starts thinning slowly through the dynamic propagation 
along with a rapid mass loss over the whole ice sheet (black 
line) mainly due to the applied perturbations. While the 
whole ice sheet starts to grow back from the 1 1th year, when 
the perturbations are switched off at a total mass loss of 
1300Gt, the interior ice sheet keeps thinning until ~300 
years due to dynamic propagation. After reaching the 
maximum mass loss of 200 Gt, the interior ice sheet also 
starts thickening. Both the interior and low-elevation regions 
reach the steady state again at ~10ka. Figure 5 illustrates 
that the strongest impact of the 1 0 year perturbations on the 
interior ice sheet occurs at ~300 years and the impact is 
sustained for lOka until full ice-sheet dynamic recovery. 
This result indicates that even if the large mass losses at the 
margins of the GIS ended today, the interior ice sheet would 
keep thinning for another 300 years and the ice sheet would 
take thousands of years to recover fully. 

CONCLUSIONS 

We used observed thinning rates below 2000 m elevation as 
the input to perturb the model to assess the inland 
propagation of the dynamic changes in the GIS. The model 
runs started from an initial present-day steady-state ice sheet 
with the perturbation applied for 10 years and then 
continued to lOOka after the perturbations were switched 
off. The results show that the dynamic propagation pene- 
trates deep into the interior ice sheet and that thinning is 
obvious within 10 years of the perturbation. The modeled 
thinning pattern over the ice sheet is in agreement with the 
observations: enhanced inland thinning mainly occurs in the 
western and eastern and, to some extent, the southeastern 


regions of the ice sheet. This implies that the strong mass loss 
since the early 2000s at the margins has had a dynamic 
impact on the entire GIS. The strongest impact for an 
isolated 1 0 year perturbation occurs at 300 years and it takes 
1 0 ka for complete dynamic recovery. We conclude that the 
interior ice response to the changes at the ice-sheet margins 
is more rapid than generally expected and the dynamic 
impact will be sustained for a long period. 
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